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ABSTRACT 

Sub- mm observations of protoplanetary disks now approach the acuity needed to measure the turbu- 
lent broadening of molecular lines. These measurements constrain disk angular momentum transport, 
and furnish evidence of the turbulent environment within which planetesimal formation takes place. 
We use local magnetohydrodynamic (MHD) simulations of the magnetorotational instability (MRI) 
to predict the distribution of turbulent velocities in low mass protoplanetary disks, as a function of 
radius and height above the mid-plane. We model both ideal MHD disks, and disks in which Ohmic 
dissipation results in a dead zone of suppressed turbulence near the mid-plane. Under ideal conditions, 
the disk mid-plane is characterized by a velocity distribution that peaks near v ~ O.lcs (where Cg is 
the local sound speed), while supersonic velocities are reached a,t z > 2>H (where H is the pressure 
scale height). Residual velocities of u w 10~^Cs persist near the mid-plane in dead zones, while the 
surface layers remain active. Anisotropic variation of the linewidth with disk inclination is modest. 
We compare our MHD results to hydrodynamic simulations in which large-scale forcing is used to 
initiate similar turbulent velocities. We show that the qualitative trend of increasing v with height, 
seen in the MHD case, persists for forced turbulence and is likely a generic property of disk turbulence. 
Percent level determinations of v at different heights within the disk, or spatially resolved observations 
that probe the inner disk containing the dead zone region, are therefore needed to test whether the 
MRI is responsible for protoplanetary disk turbulence. 

Subject headings: accretion, accretion disks — (magnetohydrodynamics:) MHD — line: profiles — 
turbulence — protoplanetary disks 
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1. INTRODUCTION 

Understanding the structure of protoplanetary disks 
is central to m odeling the phenomeno logy of Young 
Stellar Objects (iWilliams fc Ciezal [20T1I ) and to theo- 
retical studies of all phases of the formation of plane- 
tary systems. For a significant fraction of their lives, 
gas within protoplanetary disks is observed to actively 
accrete onto the central star, probably as a conse- 
quence of turbulent tran sport of angular momentum 
(jShakura fc Sunvaevlll973[ ). Globally, this accretion and 
redistribution of angular mome ntum results in evolu- 
tion o f the surface density profile (jLvnden-Bell fc Pringld 
Il974| ) , which limits the time scale for massive planet for- 
mation an d affects quantities s uch as the rate of planet 
migration (jLubow fc Idal[20Tol ). On smaller scales, tur- 
bulence sets the environ ment for planetesimal formation 
(jChiang fc Youdinll201Cl[ ) by determining both the local 
concentration and co llision velocities (fVolk et al. 1980; 
lOrmel fc Cuzzil[2007t ) of small particles that are aerody- 
namically coupled to the gas. 

Although several physical processes — including self- 
gravity and the magn etorotational instability (MRI; 
iBalbus fc Hawlevlll998f) — may in itiate disk turbulence 



(for a review, see lArmitagg 120 lit ) , existing theoretical 
studies have largely been untroubled by observational 
validation. The most widely accepted constraint on 
disk t urbulence comes from measu rements of disk life- 
times (jHaisch. Lada fc Ladal l200l|) and accretion rates 
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(jHartmann et al.l Il998[ ) , which imply that protoplane- 
tary disks around low-mass stars evolve and are dis- 
persed on Myr time scales. This observation pins down 
the angular momentum transport efficiency if the evolu- 
tion results from turbulence; the efficiency is conv e ntion - 
ally expressed in terms of a iShakura fc SunvaevI (|1973[ ) 
a K, 10^^. Generically, this level of stress within a 
fluid disk implies characteristic velocity perturbations 
y ^ Q^I'^Cs ^ 04.'"'' (w here Cs is the sound speed, e.g. 
IBalbus fc Hawlevlll998() . but this estimate is so crude as 
to be useful mainly for motivating further observations. 
Neither it, nor other constraints on a from detailed mod- 
eling of individual systems (Hucso fc Guillot 2005) pro- 
vide any information on the nature of turbulence or on 
any dependence of its properties on height above the mid- 
plane. 

Direct determination of the strength of protoplane- 
tary disk turbulence is possible by detecting the tur- 
bulent br oadening of m olecula r lines ob served in the 
infrared (iCarr. Tokunaga fc Naiital |2004[) or sub-mm 
(jHughes et al.ll2011|) . Subsonic turbulent broadening is a 
challenging quantity to measure, as protoplanetary disks 
are comprised of supersonically orbiting gas; thus, pre- 
cise measurements are needed to separate the small tur- 
bulent component from the dominant bulk rotation. Fur- 
thermore, in the inner disk, observed lines fr om the disk 
may be contaminated by outflow components (|Bast et al.l 
120111) . Nonetheless, current observations of the outer re- 
gions of disks already attain precisions comparable to 
the level {v ^ O.lc^) where a signal can plausibly be ex- 
pected. Using the Submillimeter Array {SMA) to observe 



the CO(3-2) transition. iHughes et al.l()201lD derived con- 
straints on the turbulent hnewidth in the atmosphere of 
the disks surrounding the T Tauri star TW Hya and the 
Herbig Ae star HD 16329^3- For TW Hya they placed an 
upper limit to the turbulent velocity of w < O.lCs, while 
for HD 163296 they obtained a tentative detection of tur- 
bulent broadening corresponding to w « OAcs- Although 
still preliminary, these observations provide a clear indi- 
cation that ALMA, with superior sensitivity and spatial 
resolution, will constrain disk turbulence to theoretically 
interesting levels for these and other sources. 

In this paper, our goal is to quantify the expected 
turbulent velocities in protoplanetary disks as a func- 
tion of radius and height above the mid-plane. We 
focus on low mass disks (TW Hya would be a good 
example) which are stable against self-gravity, and as- 
sume that the MRI is the sole source of turbulence. 
We compute both reference models in the ideal mag- 
netohydrodynamic (MHD) limit, and physical models in 
which the MRI is partially damped by Ohmic dissipation , 
forming a dead zone JG ammic 199J; iSano et all 120001 : 
iFromang. Terquem fc BalbusI l2002l) . Particular care is 
taken to ensure that the results are numerically con- 
verged; we use local shearing box simulations whose 
convergence w ith spatial resolution has previously been 
demo nstrated ()Davis. Stone &: Pessah|[2010t iSimon et al.l 
120111) . and we explicitly check the effect of varying the 
domain size. We also calculate purely hydrodynamic 
simulations in which turbulence is initiated through ar- 
bitrary large-scale forcing. By comparing these to the 
MHD runs, we address the question of whether the ob- 
servable properties of disk turbulence can constrain the 
underlying mechanism that initiates angular momentum 
transport. 

The plan of the paper is as follows. In §2 we describe 
the numerical simulations in ideal MHD, non-ideal MHD, 
and pure hydrodynamics that form the basis of the tur- 
bulent velocity calculation. In §3 we outline the velocity 
distribution calculation, the results of which are shown 
in §4. §5 discusses our results and their implications for 
future observations. Finally, we summarize our conclu- 
sions in §6. 

2. SIMULATIONS 

2.1. Numerical Method 

We numerically solve the equations of magnetohydro- 
dynamics (MHD) using the shearing box approximation. 
The shearing box is a model for a local, co-rotating disk 
patch whose size is small compared to the radial distance 
from the central object, Rq. This allows the construc- 
tion of a local Cartesian frame {x, y, z) that is defined in 
terms of the disk's cylindrical co-ordinates (i?, (j), z') via 
X ^ [R — Rq), y — i?o0: and z = z'. The local patch 
co-rotates with an angular velocity J7 corresponding to 
the orbital fr e quenc y at Rq, the center of the box; see 
iHawlev et al.l ()1995f ) and Fi gure [H In this fram e, the 
equations of motion become ( Hawlev et al.l [19951 ): 



-^+V-{pvv - BB)+V ip+ -B^ ] = 2qpn'^x-pn'^z-2nxpv, 



dp 
dt 



V • (pv) - 0, 



(1) 



dB 



V X {vx B) = -V X {t]V X B) 



(2) 
(3) 



where p is the mass density, pv is the momentum den- 
sity, B is the magnetic field, P is the gas pressure, and 
q is the shear parameter, defined as g = —dlnVl/ dhiR. 
We use q = 3/2, appropriate for a Keplerian disk. We 
assume an isothermal equation of state P — pc^, where 
Cs is the isothermal sound speed. From left to right, the 
source terms in equation ([2]) correspond to radial tidal 
forces (gravity and centrifugal), vertical gravity, and the 
Coriolis force. The source term in equation ([3]) is the 
effect of Ohmic resistivity, ry, on the magnetic field evo- 
lution. Note that our system of units has the magnetic 
permeability p — 1. 

Adopting this shearing box approximation allows for 
better resolution of small scales within the disk, at the 
expense of excluding global effects (th ose of scale ^ 
Ro) w hich could be physically important (|Sorathia et al.l 
120111 ). For our purposes this trade-off is worthwhile, 
because we need to numerically resolve non-ideal MHD 
terms, such as Ohmic dissipation, that play an impor- 
tant role in the structure and evolution of protoplanetary 
disks (e.g.. 'Si mon et al.ll2011[ ). 

Our simulations use Athena, a second-order accurate 
Godunov flux-conservative code for solvin g the equa- 
tions of MHD flOardiner fc Stond lMil [20081: IStone et al.l 
120081 : iStone fc Gardinedl2010( r^The numerical integra- 
tion of the shearing box equations require additions 
to the A thena algorithm, t he d etails of which can be 
found in [Stone fc Carding ([2010,^ and the Appendix of 
[Simon et al.l ( 2011t ). Brieflv. we utilize Crank-Nicholson 
differencing to conserve epicyclic motion exactly and or- 
bital advection to subtra ct off the background shear flow 
([Stone fc Gardinen 120101 ). The y boundary conditions 
are strictly p eriodic, whereas the x boundaries are shear- 
ing periodic ([Hawlev et al.[[1995l:[Simon et al.|[2011[ ). The 
vertical bou ndaries are the outfl ow boundary conditions 
described in [Simon et al.[ ([20 lit ). Finally, for simulations 
that include Ohmic resistivity, the resistive term is added 
via flrst-order in time operator splitting. We also run two 
simulations in the purely hydrodynamic limit (i.e., with 
no magnetic fields), for which we use the HLLC Rie- 
mann solver ([Toro|[1999l : [Stone et"aLl [2008() appropriate 
for non-MHD fluids. 

2.2. Runs, Parameters, and Initial Conditions 

Most of our calculations include MHD, and focus on 
the turbulent state of the MRI. The MHD simulations 
are broken down into two groups. 

The first set focuses on the ideal MHD limit, in which 
no physical dissipation is included. These simulations are 
vertically stratified, with an initial density corresponding 
to isothermal hydrostatic equilibrium. 



^ These are not "typical" sources. TW Hya is a nearby system 
with a favorable near face-on geometry, while HD 163296 has a 
very large (500 AU) disk. 



p{x,y,z) = poexp [~jp] 



(4) 



where po = 1 is the mid-plane density, and H is the scale 
height in the disk. 
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Fig. 1. — Schematic illustration of the calculation of turbulent 
velocity distributions and the relationship between a local simula- 
tion domain and a disk inclined by an angle i with respect to the 
observer. The local domain is a co-rotating patch of the larger disk, 
the size of which is small enough to approximate this domain as a 
Cartesian box. We extract the turbulent velocity from this local 
domain, appropriately averaging over time and azimuthal angle <j), 
as outlined in the text. 
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V2c 



(5) 



The isothermal sound speed, Cs = 7.07 x 10 "*, cor- 
responding to an initial value for the gas pressure of 
P„ = 5 X IQ-'^. With n = 0.001, the value for the scale 
height is iJ = 1. 

For all ideal MHD runs except for the largest domain 
calculation, the initial magneti c field confi guration is the 
twisted azimuthal flux tube of iHirose et a l. (2006), with 
minor modifications to the dimensions and the value of 
the gas to magnetic pressure ratio, /? — 2P/B^. In par- 
ticular, the initial toroidal field, By, is given by 



i?,. 




(S^ + B^^) 



HBl 



Bl^Q 



(6) 



with the toroidal /3 value Py — 100. The poloidal field 
components, B^ and B^, are calculated from the y com- 
ponent of the vector potential. 



A7, 



cos 



(ir)] 



if r < a 
if r > a 



(7) 



where r = ^/x^ + z^ and fip = 1600 is the poloidal field j3 
value. We choose a to always be one fourth of the radial 
domain size; a = Lx/4:. 

The largest domain run is initialized with a volume- 
filling toroidal field at a constant /3. We seed the MRI in 
these runs by introducing random perturbations to the 
density and velocity componentso In order to classify 
any dependence of our results on the size of the local 
region we examine, we have run several domain sizes: 
{Lx,Ly,L^) = 2H xAH X 8H, iH x 8H x 8H, 8H x 
16H X 8H, and 16H x 32H x 8H. The largest domain 

^ These init ial conditions are id entical to those of the ideal MHD 
simulations of lSimon et al.l l |2011| ). 



run has a resolution of 36 grid zones per H, and the 
resolution of each of the other runs is 32 zones per H. 

Ideal MHD is not a good approximation for most 
radii in protoplanetary disks. We have therefore run a 
second set of MHD simulations that include a height- 
dependent Ohmic resistivity ri{z), whose effect is to damp 
the MRI in regions where the resisti vity is sufficiently 
high ([Fleming. Stone fc Hawlevll2000[) . The first princi- 
ples calculation oiri{z) at different radial locations within 
the disk is difficult, because the resistivity depends on 
both the sources of ionization and on the recombination 
rate. The latter is particularly uncertain, because it is 
tied to the unkn own size distribution of small dust grains 
(|Armitagell201lD . Here, we ado pt a simple approa, c h that 
follo ws that use d previously bv lFleming fc Stoni (|2003D 
and [Turner fc S ano (200 ^. We ad opt a minimum mass 
solar nebula model (|Havashilll981 ^. and account for ion- 
ization from X-rays, cosmic rays, and the radioactive de- 
cay of ^^Al. For recombination, we consider only gas 
phase processes, and neglect dust physics. The resistiv- 
ity is related to the electron fraction Xe by. 



jy = 6.5 X 10"^a;g ^cm^s 



(|Havashilll981f) . where, assuming charge neutrality. 



Xq 



e 



FriH 



1/2 



(8) 



(9) 



Here ^ is the ionization rate, comprised of the cosmic ray 
ionization rate, 



eCR = 10 



-17 



~S„(z)/100gcm" 



-i;(,(z)/100gcm 



the X-ray ionization rate (jTurner fc Sanoir2008t ). 



(10) 



-15 



(^) 



-Sa(z)/8gcm~ 



-S6(z)/8gcm 



CxR = 2.6x10 

(11) 

and the ^^Al dec ay rate, which is constant at 4 x 10 ^^s ^ 
(|Stepinskilll992[ ). In these expressions ^aiz) and ^biz) 
are the column density lying above and below a vertical 
point z. The recombination rate, F, is 



F = 8.7 X 10"" - 



-1/2 



cm'^s 



(12) 



(jGlassgold et al.lll986( ). and nn is the number density of 
hydrogen, 



riH = 5.8 X 10 



Vau/ 



-11/4 



-^V-ff" 



(13) 



which is proportional to the gas density, p ()Wardlel l2007') . 
We have run three shearing boxes with this height- 
dependent resistivity, each of which was restarted from 
the AH X 8H x 8H ideal simulation at 100 orbits into the 
integration but with the appropriate resistivity profile 
added. The first, centered on Rq — 4AU has a large 
dead zone within ~ 2H of the mid-plane surrounded 
by two MRI active regions. The second region is cen- 
tered on Rq = lOAU, which is an intermediate region 
where the resistivity near the mid-plane is large enough 
to cause some damping of MRI turbulence but not suf- 
ficiently large to completely quench the turbulence, re- 



TABLE 1 
Shearing Box Simulations 



Label 



Domain Size 

{Lj: X Ly X Lz)H 



Resistivity? MHD/Hydro KE 



(l''Jhl/'^s)pcak 

for z > 3H 



{|''z|/Cs)pcak 

for z> ZH 



% \Vh\/c, > 1 

for z>3H 



% \V.\/Cs > 1 

for z > 3H 



Ideal-Lx2Ly4Lz8 


2x4x8 


No 


MHD 


0.03 


0.38 


0.37 


3.2 


1.9 


Ideal-Lx4Ly8Lz8 


4x8x8 


No 


MHD 


0.02 


0.54 


0.45 


5.1 


5.3 


Ideal-Lx8Lyl6Lz8 


8 X 16 X 8 


No 


MHD 


0.02 


0.61 


0.48 


9.3 


7.5 


Ideal-Lxl6Ly32Lz8 


16 X 32 X 8 


No 


MHD 


0.03 


0.66 


0.54 


13.6 


9.0 


Resistive-4AU 


4x8x8 


r]{z) at 4 AU 


MHD 


0.002 


0.56 


0.48 


9.5 


10.5 


Resistive-IOAU 


4x8x8 


rj(z) at 10 AU 


MHD 


0.02 


0.56=^, 0.54'' 


0.45=^, 0.52'' 


12. 1=^, g.i'' 


12. 7=^. 8.1*^ 


Resistive-50AU 


4x8x8 


ri{z) at 50 AU 


MHD 


0.03 


0.54 


0.48 


6.0 


7.7 


Hydro-HA 


4x8x8 


No 


Hydro 


0.07 


0.45 


0.54 


6.7 


8.3 


Hydro-LA 


4x8x8 


No 


Hydro 


0.007 


0.22 


0.18 


0.03 


0.08 



'^ Low stress state value 
^ High stress state value 

suiting in episodic bursts of mid-pla ne turbulence resem- 
bling the constant resistivity runs of lSimon et al.l ()201in . 
This dramatic variability results from the competition 
between Ohmic damping of MRI turbulence and the 
shearing of residual radial field into toroidal field of suf- 
ficient strength to reactivate the turbulence. Finally, the 
third shearing box is centered on Rq — 50AU and has 
sustained turbulence throughout the domain as the re- 
sistivity is not large enough to damp out the MRI. 
Thus, we have three different physical regimes for MRI- 



driven turbulence : one tha t resembles the classic layered 
accretion model (jGammid [1996?) . one that is relatively 
close to ideal MHD, and one intermediate regime that 
leads to large amplitude variability in turbulence lev- 
els. We should note that the radial locations of these 
regimes are subject to some uncertainty given the par- 
ticular disk model that we have adopted. If we adopted 
another model, such as a constant a disk model for ex- 
ample, then we may find the radial locations of our three 
regimes would be different. 



Finally, as one of our goals in this work is to explore 
how sensitive our derived turbulent velocity distributions 
are to the underlying mechanism for generating turbu- 
lence, we have also run forced turbulence hydrodynamic 
shearing boxes. These runs are also isothermal, verti- 
cally stratified with an initially exponential density pro- 
file (Equation |4]) , and have the same values for po, Cg, 
and rt. In these cases, we do not evolve the induction 
equation {B — 0), and we instead add a force to the 
momentum equation. 



f{x, y, z) — pA[s\a{kxx)cos(kyy)cos. 
xcos{kzz)y + sin(fc2z)2; 



\{kzz)x — cos{kxx)sm{kyy) 
(14) 



where k^ — A-k/Lx, ky — Sn/Ly^ k^ = Stt/L^, and A is 
the amplitude of the forcing. This forcing is only applied 
for \z\ < 2H. We have produced two of these calcula- 
tions, one with A — 10^'^ and one with A — 10^**. These 
calculations were performed at a resolution of 36 zones 
per H and at a domain size of 4iJ x 8H x 8H. 

Evolving the MHD simulations becomes difficult if 
there are magnetized regions of very low density, where a 
large Alfven speed results in a small timestep. Moreover, 
errors in energy make it hard to evolve regions of very 
strong field relative to gas pressure without encountering 
numerical problems. To avoid these problems, we apply 
a density fioor at a level of 10^'* of the initial mid-plane 
density throughout the physical domain in our MHD sim- 
ulations. We also include a density floor in our hydro 
simulations, which we set to 10^®. The hydrodynamic 
floor can be much lower since there is no Alfven speed 
restriction on the timestep. 



I 

Table [1] summarizes the runs, along with some basic 
properties of the turbulence that they generate. The 
ideal MHD runs are labelled with "Ideal" as a prefix and 
then the x,y, z domain size in units of H. The resistive 
runs have the prefix "Resistive" appended with the do- 
main's radial location in our model disk. Finally, the 
forced hydrodynamic runs are prefixed with "Hydro" , 
and suffixed with HA (for high-amplitude; A ~ 10^'^) 
or LA (for low-amplitude; A ~ 10^*). 

3. VELOCITY DISTRIBUTION CALCULATION 

In this work, we do not consider any radiative trans- 
fer effects or an emission model. Instead, we determine 
how the density-weighted turbulent velocity distribution 
depends on location within a protoplanetary disk and on 
the physics that we include. Although not equivalent to 
an observed turbulent line profile, the velocity distribu- 
tion gives us the probability of observing emission at a 
particular velocity shift along the line of sight. 

The line-of-sight ilos) turbulent velocity of a patch of 
disk will depend on the inclination angle of the disk i, 
and the azimuthal angle around the disk center (see 
Figure [J), 

vios = WrCos(0)sin(i) — ti0sin((/))sin(i) -I- VzCO^{i\ (15) 

where (vr^v^^Vz) is the turbulent velocity field in cylin- 
drical coordinates centered on the disk. We can rewrite 
this velocity field in terms of shearing box coordinates 
(x, y, z) as Vx = Vr, Vy = v^, and Vz = Vz\ this is a trivial 
transformation because we are interested in the magni- 
tude of the turbulent velocity fluctuations, which is the 
same in either frame. In principle, spatially resolved ob- 



servations of disks at different inclinations could yield in- 
dependent constraints on all three velocity components. 
For simplicity, we consider here just two components, a 
vertical turbulent velocity and an azimuthally averaged 
combination of v^ and Vy that corresponds to an average 
over an annulus of the disk. This "horizontal" (i.e., disk 
planar) turbulent velocity magnitude is defined as. 



p (code units) 



Vh\ 



1 

2^ 



27r 



\Vx:COs{(j)) 



i(<^)M0- (16) 



Practically, we extract v^, Vy, and \vz\ from our shearing 
box calculations, analytically average Vx and Vy as de- 
scribed above to get |w/i|, and then time-average the re- 
sulting density- weighted velocity distributions over some 
period during the saturated state. This period varies 
greatly between our various runs and was chosen based 
upon the system being in a statistically steady state and 
the horizontally averaged density being above the floor 
in the upper \z\ regions. Finally, to represent different 
line penetration depths, we calculate each distribution 
for z > 3H, z > 2H, z> H, and z> Q. 

4. RESULTS 

Before discussing the velocity distributions themselves, 
we first explore two basic diagnostics of the turbulent 
flows that are generated in the MHD and hydrodynamic 
simulations. The first diagnostic is the turbulent kinetic 
energy normalized by the gas pressure. 



KE; 



2 (P) 



(17) 



where the brackets denote a volume average (over the 
entire domain), and the overbar denotes a time-average 
over a suitable interval in which the turbulence is statis- 
tically steady. 

Table [T] displays the normalized kinetic energy for all 
simulations. The time average is done onward from orbit 
50 for all simulations except for Hydro-LA, in which it 
was done from orbit 150 onwards. With the exception 
of the strong dead zone shearing box (Resistive-4AU) , 
KE - 0.02 - 0.03 for aU MRI simulations. The time 
history of the kinetic energy in Resistive- 10 AU is highly 
variable, up to a factor of 4. Yet, the time-averaged value 
of this kinetic energy is consistent with the other MRI 
calculations. By design, the two forced hydro simulations 
bracket the MRI simulations in terms of kinetic energy. 
Hydro-HA has more kinetic energy than most of the MRI 
simulations by a factor of ~ 3, whereas Hydro-LA has 
less kinetic energy by about the same factor. 

While we have established that we can force turbu- 
lence to roughly the same kinetic energy amplitude as the 
MRI-driven cases, it is worth asking what the structure of 
the forced turbulence is. Do the forced turbulence runs 
resemble their forcing functions at late times, or does 
a different structure emerge? In Figure [2l we plot the 
density in the mid-plane of the Resistive-50AU run and 
the Hydro-HA run at 100 orbits into the evolution. The 
two runs are visually nearly indistinguishable, and even 
in the forced turbulence case, there exist spiral density 
waves that propagate through the domain. Indeed, the 
auto-correlation function for the gas density (|Guan et al.l 
[2009. : .Nelson fc Grcsscli2010, ) returns a density structure 




Fig. 2. — Snapshot of the mid-plane gas density at 100 orbits. 
The left panel is from run Resistive-50AU and represents an MRI 
calculation. The right panel is from Hydro-HA and is a forced tur- 
bulence simulation. Both calculations show the presence of spiral 
density waves. 

that is very similar between the two runs. This simi- 
larity may be a result of the choice of forcing function 
for the hydro calculations; if we had chosen some other 
forcing function, perhaps these density waves would not 
exist or would look different t o the MRI case. How- 
ever, [Hememan^&iTaEiloiio^ (|2009) suggest that these 
waves can be generally produced by disk turbulence, not 
necessarily restricted to that which is MRI-driven. In 
this respect, our hydro calculations are a representation 
of potential forms of disk turbulence that produce these 
waves other than the MRI. 

Figure [3] displays the density-weighted turbulent veloc- 
ity distribution for several integration depths and shear- 
ing box domain sizes, all in the ideal MHD limit. The 
dashed lines are \vz\ and the sohd lines are \vh\. The 
most striking feature of these plots is the rapid increase 
in the velocity of the peak of the distribution as one 
moves higher in the disk. For z > (upper half of the 
disk) the distribution peaks at about 10% of the sound 
speed, but this velocity increases to about 50% of Cg for 
z > ?>H . The width of the distributions is quite large; for 
the AH X %H x 8H domain, there is a ~ 90% probability 
that |wh|/cs lies between 0.1 and 1. There does not ap- 
pear to be a strong difference between the \vz\ and \vh\ 
distributions, with the latter being slightly more sharply 
peaked and at a slightly higher jt'l/cs than the former. 
This suggests that the inclination of the disk will only 
weakly play a role in the observed turbulent velocity. Fi- 
nally, convergence of the peak velocity with domain size 
appears to have been attained for the AH x 8H x 8H 
domain; this suggests that all the essential physics in- 
volved in setting the magnitude of velocity fiuctuations 
is captured by this intermediate-sized domain. 

Another interesting feature to note is the non- 
negligible supersonic velocity component to the distri- 
bution above 3-ff. Integrating over the distribution for 
|w|/cs > 1 yields roughly 10% of the turbulent velocity 
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Fig. 3. — Turbulent velocity distributions in the ideal MHD 
calculations. Each panel corresponds to a different local domain 
size, and the different colors in each figure correspond to different 
depths over which the distribution is calculated, as labeled. The 
dashed lines are the vertical turbulent velocity lii^l/cs, and the 
solid lines are the azimuthally averaged disk planar velocity |i)(,|/cs. 
Most of the turbulent velocities are in the range jfl/cs ~ 0.1 — 1. 



being supersonic for box sizes larger than 8i/ x 16iJ x 8i/. 
The smaller domains have smaller supersonic compo- 
nents: -- 5% and ~ 1% for the AH x 8H x 8H and 
2H X AH X 8H, respectively. The origin of these super- 
sonic velocities is presumably the steepening of initially 
subsonic waves as the density decreases away from the 
mid-plane. Similar physical effects have b een seen in 
many prior simula.tions of stratified disks ([Stone et al.l 
[l99a iFlock et al.l [20I1I : IBeckwith et all 120 111) . Along 
with the recently studied dissipation of current sheets in 
disk corona (Hirosc & Turner 2011), shock heating from 
these supersonic motions could potentially play an im- 
portant role in the thermodynamic properties of disk at- 
mospheres. We will further explore these thermodynamic 
issues in future publications. We include the peak of the 
|w|/cs distribution and the percentage of the distribution 
with supersonic velocities for z > 3H in Table [TJ 

In Figure [4] we plot the same velocity distributions for 
non-ideal (resistive) shearing boxes computed at differ- 
ent radial locations. Since the simulation conducted at 
10 AU is highly variable, we show results that corre- 
spond both to the high stress turbulent state, and to 
the low stress state. Considering first the shearing box 
centered on 4 AU, it appears that as one probes regions 
closer to the mid-plane, the turbulent velocity fiuctu- 
ations drop dramatically, with a peak in the distribu- 
tion at \v\/cs ~ 0.02. This is not surprising since the 
dead zone in this run extends to about ±2H, and the 
velocity fiuctuations induced by the active layers a.ppear 
within the dead zo ne region fe.g.. lFleming fc Stonell2003l : 
iSimon et al.ll201l]) . Above 3H the velocity distribution 
is very similar to the ideal MHD cases, including the 
presence of supersonic velocities. 

Moving outward to near the outer edge of the dead 
zone, the 10 AU simulation yields a velocity distribution 
that varies slightly, depending on whether or not the sys- 
tem is in the "high state" or the "low state". The low 
state appears to be intermediate in the velocity distri- 
bution between the 4 AU and ideal MHD cases, whereas 
the high state resembles the ideal MHD distribution more 
closely. In both states, the velocity distribution near the 
disk surface peaks around |w|/cs ~ 0.5 with a substan- 
tial supersonic tail, again agreeing with the other simula- 
tions. Finally, the shearing box centered on 50 AU has a 
distribution very similar to the ideal MHD case, consis- 
tent with the notion that the resistivity is small enough 
at this radius to not significantly affect the MRI. 

Figure [5] displays the velocity distributions for the two 
forced hydro runs. Hydro-HA has a distribution quite 
similar to those in the MRI calculations. However, there 
is a weaker dependence of the turbulent velocity on the 
height from the mid-plane; the peak |w|/cs values lie be- 
tween 0.2 and 0.5. Again, the peaks of the |wz|/cs and 
jw/il/cg distributions are quite similar. There is also a 
significant supersonic component to the z > SH distri- 
bution; ~ 7-8% of the distribution has \v\/cs > 1. The 
distribution peaks for Hydro-LA are lower, which is not 
surprising since there is less kinetic energy in this run. 
However, despite an order of magnitude difference in the 
saturated kinetic energies, the peak velocity for z > 3H 
is only a factor of 2.5 lower in Hydro-LA than in Hydro- 
HA. The mid-plane velocities are significantly lower in 
Hydro-LA than in Hydro-HA, however. These results 
suggest that even when forced with a lower amplitude. 
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Fig. 4. — Turbulent velocity distributions for the resistive MHD 
calculations. The lower two panels are the box at 10 AU in a 
state characterized by low turbulent stress (top) and high turbu- 
lent stress (bottom) . The different colors in each figure correspond 
to different depths over which the distribution is calculated, as la- 
beled. The dashed lines are the vertical turbulent velocity |dz|/cs, 
and the solid lines are the azimuthally averaged disk planar veloc- 
ity \vfi\/cB. The 4 AU calculation shows the presence of a strong 
dead zone, as the velocity distribution peaks at a much lower value 
towards the mid-plane. 
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Fig. 5. — Turbulent velocity distributions for the forced hydro 
calculations. The top panel is the run forced with amplitude 10^'' 
and the bottom panel has forcing with amplitude 10~ . The differ- 
ent colors in each figure correspond to different depths over which 
the distribution is calculated, as labeled. The dashed lines are 
the vertical turbulent velocity liizl/cs, and the solid lines are the 
azimuthally averaged disk planar velocity |w(,|/cg. 

the turbulent velocities can steepen significantly in the 
lower density regions away from the mid-plane. Taken 
together, the characteristic velocities in the forced hydro 
runs are not very different than those in the MRI cases. 

5. DISCUSSION AND IMPLICATIONS FOR OBSERVATIONS 

This paper represents a step toward making a direct 
connection between the simulated properties of turbu- 
lent protoplanetary disks and actual observations of these 
systems. To this end, we have presented a series of calcu- 
lations with varying physics from which we extracted the 
turbulent velocity distribution. The simulations do not 
explicitly predict actual observables, as that would re- 
quire the inclusion of radiation physics in one form or an- 
other. However, our results do have several implications 
for the nature of disk turbulence in low mass protoplan- 
etary disks that could potentially be tested with future 
observations, particularly those made with ALMA. 

The first implication is that if turbulence is driven 
solely by the MRI, the turbulent linewidth ought to vary 
strongly as a function of both radius and height above 
the mid-plane. In regions of the disk where non-ideal ef- 
fects are small, we predict mid-plane velocities that peak 
near 0.1 times the local sound speed. The characteristic 
turbulent velocities increase with height, such that sig- 
nificant regions of transonic flow occur in the atmosphere 
at z > 3H. This prediction is in at least rough agree- 
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Fig. 6. — Time- and horizonally-averaged vertical density profiles 
for a subset of the shearing box simulations. The time average was 
done from orbit 50 onward. The density is in code units (pz=o ~ !)> 
and 2 is in units of H . The green curves correspond to the ideal 
MHD simulations, and the solid green curve is Ideal-Lx2Ly4Lz8, 
the dashed green curve is Ideal-Lx4Ly4Lz8, the dotted green curve 
is Ideal-Lx8Ly8Lzl6, and the triple-dot dash green curve is Ideal- 
Lxl6Ly32Lz8. The blue line corresponds to Resistive-4AU, the 
black line is Resistive-50AU, and the red line is Hydro-HA. The 
forced hydro run has a nearly Gaussian density profile, whereas the 
MHD calculations have a Gaussian profile for \z\ < 2H, outside of 
which the density gradient flattens out. 
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Fig. 7. — Peak of the turbulent velocity distribution versus frac- 
tion of the total surface density. The squares are the planar velocity 
distribution peaks, and the asterisks are the vertical velocity peaks. 
The black symbols are the peak velocities for the MRI calculation 
at 50 AU (Resistive-50AU), blue is the MRI calculation at 4 AU 
(Resistive-4AU), and red is the higher amplitude forced hydrody- 
namic turbulence run (Hydro-HA) . As one probes vertically deeper 
into the disk, the peak velocity decreases. 

ment with the observational meas urement of tu r bulent 
broadening in HD 163296, where iHughes et all ()2011f ) 
infer turbulent line widths consistent with a few tenths 
of the sound speed. The observational results for this 
system are then consistent with MRI-driven turbulence 
actually occurring in the outer disk. The upper limit for 
TW Hya, on the other hand, is at best marginally con- 
sistent with the MRI prediction for near-ideal conditions 
- pushing that limit lower has the potential to provide a 



stringent test of our models. 

If we assume, on theoretical grounds, that the MRI is 
the only viable source of turbulence in low mass disks, 
then the agreement between our simulations and the ob- 
servational results for HD 163296 is mildly encouraging. 
That is, MRI-driven turbulence produces the "right" an- 
swer. This order of magnitude level of agreement, how- 
ever, is nowhere near discriminating enough to exclude 
the possibility that some (unspecified) hydrodynamic in- 
stability is responsible for the observed turbulent broad- 
ening. We have run purely hydrodynamic disk models 
that are driven to a turbulent state via large scale forc- 
ing, and these disks also show velocity fluctuations of the 
order of tenths of the sound speed, along with a velocity 
gradient between the mid-plane and the surface, where 
transonic velocities can be produced. 

This raises the question, can observations distin- 
guish between purely hydrodynamic turbulence and that 
driven by the MRI? Of course, one option is to ob- 
serve the magnetic field structure and amplitude (or lack 
thereof) in the disk i tself, but this is cu rrently exceed- 
ingly difficult (but see iHughes et al.l[2009l ) . If one cannot 
observe the field directly, the next best option is to ob- 
serve a secondary effect of the magnetic field such as the 
turbulent velocity. However, as noted already, a single 
measurement of the turbulent velocity is certainly insuf- 
ficient to distinguish between MHD and hydrodynamic 
drivers of turbulence. An arbitrary forcing of the fluid 
can still yield turbulent velocities that are more or less 
indistinguishable from those produced via the MRI. 

The prospects for distinguishing between different 
sources of turbulence are better if observations are able to 
probe either different heights within the disk (by exploit- 
ing multiple molecular species), different radii, or both. 
In a real disk, of course, there is no immediately accessi- 
ble observable that isolates the turbulence at a particular 
physical height above the mid-plane. Rather, the observ- 
able is the degree of broadening of a given molecular line 
produced in a region of the disk where the temperature, 
density and chemistry yield sufficient emissivity, and the 
optical depth is not too high. Thus, where in z a par- 
ticular line is emitted will depend on the vertical struc- 
ture of the disk itself. We find that this structure differs 
significantly depending upon whether the disk is magne- 
tized or not. Consider Fig. [6l which shows the time- and 
horizontally-averaged vertical density profile. The red 
curve is one of the forced hydro turbulence cases, and 
the other curves are a subset of the MRI-driven turbu- 
lence runs. There is an obvious and significant difference 
between the density profiles at large \z\. The density de- 
parts significantly from Gaussian in the MRI cases. The 
reason for this is that for \z\ > 2_ff, magnetic pressure 
dominates over gas pressure, and the gradient in mag- 
netic pressure helps to support the gas against gravity. 
Thus, the gas pressure and density have a shallower slope 
in these regions. This magnetic and gas pressure struc- 
ture i s consistent with p revious shearing box simulations 
(e.g., iHirose et al1l2006( ). 

The difference in vertical structure between a magne- 
tized and non-magnetized disk results in a distinct differ- 
ence in the the profile of turbulent velocity with column 
density. This is shown in Figure [71 which plots how the 
characteristic turbulent velocity changes as a function of 
fractional column, defined as 
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where p{z) is the time- and horizontaUy averaged gas 
density (the time average is done from orbit 50 onwards). 
There are two features to note in this plot. The first is 
that, in general, the turbulent velocity decreases as one 
probes deeper into the disk. This result was discussed 
above in the context of the velocity distributions, and is 
simply the result of velocity steepening in lower density 
regions. The second feature is that the v/cs ^ 0.5 values 
obtained in the upper disk regions can be found at a 
lower AS/E (by about an order of magnitude) in the 
hydro case versus the MHD cases. This suggests that 
the column depth to which a particular line can probe 
may be very useful in determining the density structure 
away from the disk mid-plane and thus can constrain the 
turbulence mechanisms. 

Furthermore, the radial dependence of the velocity gra- 
dient with distance from the mid-plane may also be use- 
ful. While the general trend of decreasing |f |/Cs with 
height is robust in all of our calculations, the presence 
of the MRI dead zone dramatically changes this gradi- 
ent. In particular, as Figure [7] shows (blue points), the 
presence of a magnetically dead zone is quite obvious as 
the turbulent velocities drop well below O.lcs near the 
mid-plane region, but the active layers above and below 
the mid-plane, in combination with steepening, produce 
|w|/cs ~ 0.5. Thus, if one were to probe different depths 
into the disk and find a dramatic decrease in turbulent 
velocity towards the mid-plane, this would be strongly 
indicative of a dead zone region. 

6. CONCLUSIONS AND UNCERTAINTIES 

Our conclusions about the turbulent properties of low 
mass protoplanetary disks are: 

1. Characteristic turbulent velocities are ^ (O.l-l)cs 
for fully turbulent regions, in rou gh agreement 
with observations made to date fe.g.. lHughes et all 
120111 ). These characteristic velocities are reason- 
ably robust to variations in numerical (changes 
to local domain size) and physical (locations in a 
model disk) parameters. 

2. Turbulent velocity increases away from the disk 
mid-plane due to steepening. In the upper region 
of the disk (|2;| > 3H), the velocity distribution 
peaks around 0.5cs and has a significant (^ 10%) 
supersonic component. As one probes towards the 
mid-plane, |i'|/cs ^ 0.1 is typical of fully turbulent 
disks. 

3. In calculations with an MRI dead zone near the 
mid-plane, the characteristic turbulent velocities 
are ~ O.OlCg within the dead zone. In principle, 
with an improvement in sensitivity, observations 
that probe different depths could see the presence 
of the dead zone as velocities drop from ^ O.l-lcg 
to ~ O.OlCs. 

4. The density structure for \z\ > 2H is signifi- 
cantly different in the MRI versus purely hydro- 



dynamic cases, which could have potential impli- 
cations for the observed turbulent linewidths if dif- 
ferent depths can be probed. 

5. The vertical and planar velocity distributions are 
quite similar, suggesting that turbulent linewidths 
will only weakly be dependent on the inclination 
angle. 

Our predictions for the distribution of turbulent veloc- 
ity in MRI-active disks suffer from a number of uncertain- 
ties. First, we have chosen to only focus on one non-ideal 
MHD effect, namely Ohmic resistivity. Other non-ideal 
effects - ambipolar diffusion and t he Hall term - are also 
important in protoplanetary disks ('Kunz fc Balbusll2004l : 
IBai fc Stonel[20Tll : IWardlefc Salmeron 20111). Ambipo- 
lar diffusion, in particular, is important in low density 
regions, and may affect the properties of turbulence in 
the most observationally accessible location - the disk 
atmosphere at large radius. Moreover, the resistivity 
that we have employed neglects dust physics. We also 
caution that some of the most striking qualitative trends 
that we observe are linked to the steepening of waves 
near the disk surface. Wave propagation in disks is 
known to depend upon the vertical thermal structure 
(|Bate et al.ll2002[ ). and hence the isothermal structure 
that we have assumed may not always be adequate. Even 
before adding a treatment of the radiation physics, these 
limitations imply that there remains much work to be 
done to further constrain the turbulent velocities in sim- 
ulations of MRI-active disks. 

While our results coupled with recent observations pro- 
vide mild support for the model of MRI-driven angular 
momentum transport, the calculations that we have pre- 
sented here have not been able to identify a strong ob- 
servational discriminant between MRI-driven and purely 
hydrodynamic turbulence. Quite precise measurements 
of the turbulence as a function of height will be needed to 
tell one from the other on purely observational grounds. 
In fact, if sufficiently high spatial resolution observations 
of the inner disk regions reveal the presence of a dead 
zone region, then this would present very strong support 
for the MRI driving disk turbulence. 

Finally, we reiterate that we have considered only ar- 
bitrary hydrodynamic forcing, rather than setting up 
known physical drivers of turbule nce, such as self-gravity 
or even convection (L esur fc Ogi lvic 2010). By design, 
the average kinetic energy in the purely hydrodynamic 
simulations nearly equals that of the MRI simulations. 
If it could be established, theoretically, that hydrody- 
namic drivers of turbulence were necessarily weaker than 
the MRI, a single measurement of the turbulent velocity 
would then distinguish between the two. If, on the other 
hand, we treat the strength of hydrodynamic turbulence 
as a free parameter, then our results suggest that the 
vertical variation of turbulent velocities in the hydrody- 
namic and MHD limits can have a qualitatively similar 
trend. Of course, the physical mechanisms that might 
initiate hydrodynamic turbulence without arbitrary forc- 
ing could, in principle, imprint distinctive characteristics 
into the observable turbulent velocity field, which would 
allow them to be distinguished from the MRI more read- 
ily. To test this, it would be useful to repeat the analysis 
presented here for disks in which these other sources of 
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turbulence are active. These calculations are currently 
underway and will be presented in future work. 
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